      SUBROUTINE DER(T,X,DX)
C  VERSION OF 4/18/28
C  PURPOSE   
C    COMPUTE DERIVATIVES OF THE EQUATIONS OF MOTION
C  INPUT ARGUMENTS
C    T      = CURRENT JULIAN TIME (SEC)
C    X      = 6-D CURRENT CARTESIAN STATE, X,Y,Z,XD,YD,ZD, (KM,KM/SEC)
C  INPUT COMMON BLOCKS
C    OPTION, TIME, PLTCON, ATMCON, SPCCON, SUNCON, MUNCON, HARMON
C  OUTPUT
C    DX     = 6-D DERIVATIVES OF POSTION AND VELOCITY  
C  CALL SUBROUTINES  
C    LEGEND, ANGLES, EPHEM, XTHIRD, SETTHD, DENS76
C  REFERENCES
C    JPL EM 312/87-153, 20 APRIL 1987
C  ANALYSIS  
C    JOHNNY H. KWOK - JPL  
C  PROGRAMMER
C    JOHNNY H. KWOK - JPL  
C  PROGRAM MODIFICATIONS 
C    NONE  
C  COMMENTS  
C    THE HARMONICS ARE DIMENSIONED FOR 40X40, WHERE COEFS. ARE 
C    STORED AS CIJ=C(I+1,J+1), ETC. FOR HIGHER FIELDS, CHANGE  
C    DIMENSION ON P, CN, SN, TN, C, S. 
C   
C  MODIFICATION HISTORY
C   Changed TIME common block name to TIMECMN, to eliminate linker 
C   linker warnings.  B. A. Franz, GSC, November 14, 1997.
C   


      IMPLICIT DOUBLE PRECISION (A-H,O-Z)   
      DIMENSION X(6),DX(6)  
      DIMENSION P(41,41),CN(41),SN(41),TN(41)   
      DIMENSION XS(3),XT(3),XM(3),XN(3),VB(3)
      COMMON/OPTION/L,M,IRES,ISUN,IMOON,IEPHEM,IDRAG,IDENS,ISRP,IORB
     1  ,IPRINT,INODE,IPLOT
      COMMON/TIMECMN/TI,TF,TR
      COMMON/PLTCON/GE,RE,RATE,PM,AJ2,ELLIP,RATM
      COMMON/ATMCON/RDENS,RHT,SHT,ALTMAX,WT
      COMMON/SPCCON/AREAD,AREAS,SCMASS,CDRAG,CSRP
      COMMON/SUNCON/GS,ES(7),ET(7)
      COMMON/MUNCON/GM,EM(7),EN(7)
      COMMON/HARMON/C(41,41),S(41,41)   
      DATA TPI/6.283185307179586D0/
      DATA ZERO,HALF,ONE,TWO,TRHALF/0.D0,0.5D0,1.D0,2.D0,1.5D0/
      R2=X(1)**2+X(2)**2+X(3)**2
      R1=DSQRT(R2)  
      R3=R1*R2
      DX(1)=X(4)
      DX(2)=X(5)
      DX(3)=X(6)
      DO 10 I=4,6   
   10 DX(I)=ZERO
      IF (L.EQ.0.AND.M.EQ.0) GO TO 200
      IF (L.EQ.2.AND.M.EQ.0) GO TO 100
C
C  COMPUTE DERIVATIVE DUE TO SPHERICAL HARMONICS OF DEGREE L AND M
C
      SR2=X(1)**2+X(2)**2   
      SR1=DSQRT(SR2)
      SPHI=X(3)/R1  
      PHI=DASIN(SPHI)   
      AMDA=DATAN2(X(2),X(1))-DMOD(PM+RATE*(T-TR),TPI)   
      IM=M  
      IF (M.LT.L) IM=M+1
      CALL LEGEND(L,IM,SPHI,P)  
      CALL ANGLES(M,AMDA,PHI,CN,SN,TN)
      E1=ZERO
      E2=ZERO
      E3=ZERO
      D1=RE/R1  
      D2=D1
      DO 110 IL=2,L 
      F1=ZERO
      F2=ZERO
      F3=ZERO
      IL1=IL+1  
      DO 120 IM=0,IL,IRES
      IF (IM.GT.M) GO TO 130
      IM2=IM+2
      IM1=IM+1
      D3=C(IL1,IM1)*CN(IM1)+S(IL1,IM1)*SN(IM1)  
      F1=F1+D3*P(IL1,IM1)   
      IF (IM2.LE.IL1) F2=F2+D3*(P(IL1,IM2)-TN(IM1)*P(IL1,IM1))  
      IF (IM2.GT.IL1) F2=F2-D3*TN(IM1)*P(IL1,IM1)   
      IF (IM.NE.0)  
     XF3=F3+IM*(S(IL1,IM1)*CN(IM1)-C(IL1,IM1)*SN(IM1))*P(IL1,IM1)
  120 CONTINUE  
  130 CONTINUE  
      D2=D2*D1  
      E1=E1+D2*IL1*F1   
      E2=E2+D2*F2
      E3=E3+D2*F3
  110 CONTINUE  
      D3=GE/R1
      E1=-E1*D3/R1  
      E2=E2*D3  
      E3=E3*D3  
      D1=E1/R1  
      D2=X(3)/(R2*SR1)*E2   
      D3=E3/SR2 
      DX(4)=DX(4)+(D1-D2)*X(1)-D3*X(2)  
      DX(5)=DX(5)+(D1-D2)*X(2)+D3*X(1)  
      DX(6)=DX(6)+D1*X(3)+SR1*E2/R2
      GO TO 200
C
C  COMPUTE DERIVATIVE DUE TO J2 ONLY.  IF ONLY J2 IS NEEDED, THIS IS A
C  MUCH SIMPLER SET OF EQUATIONS TO USE.  NOTE THAT AJ2 IS -C20
C
  100 CONTINUE
      R5=R2*R3
      D1=-TRHALF*AJ2*RE*RE*GE/R5
      D2=ONE-5.D0*X(3)**2/R2
      DX(4)=DX(4)+X(1)*D1*D2
      DX(5)=DX(5)+X(2)*D1*D2
      DX(6)=DX(6)+X(3)*D1*(D2+TWO)
  200 CONTINUE
C
C *** SETUP LUNI-SOLAR POSITION IF IEPHEM.EQ.1
C
      IF (IEPHEM.EQ.0) THEN
        TRF=TR
      ELSE
        IF (ISUN.EQ.1.OR.IMOON.EQ.1.OR.ISRP.EQ.1)
     1  CALL EPHEM(ISUN,IMOON,ISRP,T,ES,EM)
        TRF=T
      ENDIF
C
C  COMPUTE DERIVATIVE DUE TO THE SUN
C
      IF (ISUN.EQ.0) GO TO 300  
      CALL XTHIRD(T,TRF,ES,ET,XS) 
      RS=DSQRT(XS(1)**2+XS(2)**2+XS(3)**2)
      RS3=RS**3
      DO 210 I=1,3
  210 XT(I)=X(I)-XS(I)
      RT=DSQRT(XT(1)**2+XT(2)**2+XT(3)**2)
      RT3=RT**3
      DO 220 I=1,3
  220 DX(I+3)=DX(I+3)-GS*(XT(I)/RT3+XS(I)/RS3)
  300 CONTINUE
C
C  COMPUTE DERIVATIVE DUE TO SOLAR RADIATION PRESSURE
C
      IF (ISRP.EQ.0) GO TO 400
      IF (ISUN.EQ.0) THEN
        CALL XTHIRD(T,TRF,ES,ET,XS)
        RS=DSQRT(XS(1)**2+XS(2)**2+XS(3)**2)
      ENDIF
      DO 310 I=1,3
  310 XS(I)=XS(I)/RS
C
C  CHECK FOR SHADOW CONDITION
C
      RDRS=X(1)*XS(1)+X(2)*XS(2)+X(3)*XS(3)
      IF (RDRS.GE.ZERO) GO TO 320
      THETA=DACOS(RDRS/R1)
      IF (R1*DSIN(THETA).LT.RATM) GO TO 400
  320 CONTINUE
      DO 330 I=1,3
  330 DX(I+3)=DX(I+3)-CSRP*AREAS*XS(I)/SCMASS
C
C  COMPUTE DERIVATIVE DUE TO THE MOON
C
  400 CONTINUE
      IF (IMOON.EQ.0) GO TO 500
      IF (IEPHEM.EQ.1) CALL SETTHD(EM,EN)
      CALL XTHIRD(T,TRF,EM,EN,XM) 
      RM=(XM(1)**2+XM(2)**2+XM(3)**2)**TRHALF
      DO 410 I=1,3
  410 XN(I)=X(I)-XM(I)  
      RN=(XN(1)**2+XN(2)**2+XN(3)**2)**TRHALF
      DO 420 I=1,3  
  420 DX(I+3)=DX(I+3)-GM*(XN(I)/RN+XM(I)/RM)
C
C  COMPUTE DERIVATIVE DUE TO DRAG
C 
  500 CONTINUE
      IF (IDRAG.EQ.0) GO TO 600
      IF (ELLIP.EQ.ZERO) THEN
        EALT=R1-RE
      ELSE
        PHI=DASIN(X(3)/R1)
        EALT=R1-DSQRT(RE**2*(ONE-ELLIP**2)/(ONE-(ELLIP*DCOS(PHI))**2))
      ENDIF
      IF (EALT.GT.ALTMAX) GO TO 600
C
C  COMPUTE DENSITY AT ALTITUDE.  THIS PORTION CAN BE REPLACED BY MORE
C  SOPHISCATED MODEL
C
      IF (IDENS.EQ.0) DENS=RDENS*DEXP((RHT-EALT)/SHT)
      IF (IDENS.EQ.1) CALL DENS76(EALT,DENS)
      DENS=DENS*WT
C
      VB(1)=X(4)+X(2)*RATE
      VB(2)=X(5)-X(1)*RATE
      VB(3)=X(6)
      VBM2=VB(1)**2+VB(2)**2+VB(3)**2
      VBM1=DSQRT(VBM2)
      DO 510 I=1,3
  510 VB(I)=VB(I)/VBM1
      DO 520 I=1,3
  520 DX(I+3)=DX(I+3)-HALF*CDRAG*AREAD/SCMASS*DENS*VBM2*VB(I)
C
C  ADD THE TWO BODY DERIVATIVE TOO
C
  600 CONTINUE
      DO 610 I=1,3  
  610 DX(I+3)=DX(I+3)-GE*X(I)/R3
      RETURN
      END   
